Model and Simulations

This is an analysis of the PTSD model “complete” simulations (“Simulation 3”).

Trajectories to simulate

One of the goals of the model is to simulate therecovery trkjaectories in PTSD. The classic distinction is between four different trajectories, Resilient, Recovery, Chronic, and Delayed. These trajectories can be stylized as such

time <- -20:60

resilient <- rep(0, length(time))
chronic <- tanh(time/2)
delayed <- time**2 / max(time**2)
#recovery <- chronic - delayed
recovery <- chronic + ((60 - time)**2 / max(time**2)-1)
#recovery <- recovery/max(recovery)

resilient[time <= 0] <- 0
chronic[time <= 0] <- 0
delayed[time <= 0] <- 0
recovery[time <= 0] <- 0

wideal <- tibble(Day = time, 
                Chronic = chronic,
                Delayed = delayed,
                Recovery = recovery,
                Resilient = resilient)
        
ideal <- wideal %>%
  pivot_longer(cols=c("Chronic", "Delayed", "Resilient", "Recovery"),
               values_to = "CTraumatic",
               names_to = "Trajectory")

ideal$Trajectory <- factor(ideal$Trajectory,
                           levels = c("Recovery", "Delayed", "Resilient",  "Chronic"))

ggplot(ideal, 
       aes(x=Day, 
           y=CTraumatic, 
           col=Trajectory, 
           fill=Trajectory)) +
  #geom_smooth(method="loess", span=0.2) +
  geom_line(size=2, alpha=0.75) +
  ggtitle("Theoretical Recovery Trajectories") +
  scale_color_d3() +
  ylab("Number of Memory Intrusions") +
  geom_vline(data=tibble(), 
             aes(xintercept=-0.35)) +
  scale_fill_d3() +
  theme_pander()  +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank())

Underlying memory model

The underlying memory model is John Anderson’s ACT-R.

# Simulate ACT-R

decay <- function(time, t0=0, rate=.5) {
    if (time <= t0) {
      NA
    } else {
      (time - t0)**(-rate)
    }
}

vdecay <- Vectorize(decay)

time <- seq(1, 100, 1)

t1 = vdecay(time, t0=1)
t2 = vdecay(time, t0=20)
t3 = vdecay(time, t0=55)
t4 = vdecay(time, t0=65)


df <-data.frame(Time=time, 
           trace1 = vdecay(time, t0=1),
           trace2 = vdecay(time, t0=20),
           trace3 = vdecay(time, t0=55),
           trace4 = vdecay(time, t0=65))


ldf <- df %>%
  # mutate(Total = if_else(is.na(trace1), 0, trace1) +
  #          if_else(is.na(trace2), 0, trace2) +
  #          if_else(is.na(trace3), 0, trace3) +
  #          if_else(is.na(trace4), 0, trace4))  %>%
  pivot_longer(cols = c("trace1", "trace2", "trace3", "trace4"),
               names_to = "Trace", values_to="Activation")

ldf$Trace <- fct_recode(ldf$Trace, 
                        "Trace 1" = "trace1",
                        "Trace 2" = "trace2",
                        "Trace 3" = "trace3",
                        "Trace 4" = "trace4")

ldf <- ldf %>%
  add_column(Source="Individual Traces")

t1[is.na(t1)] <-0
t2[is.na(t2)] <-0
t3[is.na(t3)] <-0
t4[is.na(t4)] <-0
total <- log(t1+t2+t3+t4)

totaldf <- tibble(Time=time, 
                  Activation = total,
                  Trace = "Memory",
                  Source="Memory")

ldf <- rbind(ldf, totaldf)

ldf$Source <- factor(ldf$Source, levels = c("Memory", "Individual Traces"))

ggplot(ldf, 
       aes(x=Time, y=Activation, col=Trace)) +
  geom_line(size=1,
            alpha=1) +
            #linetype="dashed") +
  # stat_summary(geom="line", fun = sum, 
  #              col="black", 
  #              alpha=0.5,
  #              position = position_dodge()) +
  # #scale_y_log10() +
  facet_wrap(~ Source, ncol=1,
             scales = "free_y") +
  #ylim(0, 1.6) +
  #ylab(expression(paste("Availability of Memory ", italic(m))))+
  ylab("Activation") +
  scale_color_viridis("", option="plasma", discrete=T, end = .8, direction=-1) +
  theme_pander()

## Warning: Removed 141 row(s) containing missing values (geom_path).

Neuroiological interpretation

brain_model <- tibble(
  region = c("parahippocampal", "entorhinal", "medial orbitofrontal",
             "superior frontal", "paracentral", "precuneus"),
  component = c("LTM", "Emotional Intensity", "Retrieval Control",
                "Context", "Context", "Context")
)

#brain_model$component <- as_factor(brain_model$component)

brain_model %>%
#  group_by(component) %>%
  ggplot(aes(fill = component)) +
  geom_brain(atlas = dk, 
             col="white",
             show.legend = T
             ) +
  labs(fill="Model\nComponent") +
  #scale_fill_lancet(na.value="lightgrey") +
  scale_fill_manual(na.value="lightgrey", 
                     values=c("dodgerblue2", "goldenrod1", "firebrick3", 
                              "lightskyblue3", "lightskyblue3", "lightskyblue3"),
                     breaks = brain_model$component) +
  coord_sf(xlim = c(750, 1050)) +
  theme_brain2(text.family = "sans",
               text.colour="black")

## merging atlas and data by 'region'

Event Distribution

The model lives in a sumulated world in which memories are automatically retrived in response to changes in the environment. These changes, called events occur probabilitistically, with Gamma distrib a predermined frequen

The distribution of Gamma events is this:

Frequency: 2.1, 175 (shape, scale)
Rumination, 6, 100

t <- 0:(24*60)
E <- dgamma(t, shape=2.1, rate=1/175)
R <- dgamma(t, shape=6, rate=1/100)

dists <- data.frame(Time=t, Events = E, Rumination = R) %>%
  pivot_longer(names_to = "Type", 
               values_to = "Probability",
               cols=c("Events", "Rumination"))

ggplot(dists, aes(x=Time, y=Probability, col=Type, fill=Type)) +
  geom_line() +
  geom_ribbon(data=dists, aes(x=Time, ymax=Probability, ymin=0), alpha=0.25) +
  scale_x_continuous(breaks=60*c(0, 4, 8, 12, 16, 20),
                     limits=c(0, 24*60),
                     labels=c("8:00am", 
                              "12:00pm",
                              "4:00pm",
                              "8:00pm",
                              "12:00am",
                              "4:00am")) +
  scale_color_aaas() +
  scale_fill_aaas() +
  ggtitle("Probability of Retrievals During the Day") +
  theme_pander()

Load and Transform the Data

First, we are going to load the “aggregated” version of the simulations, with the events already aggregated by day.

CREATE_AGGREGAGATED = F

if (CREATE_AGGREGAGATED) {
  d <- read_csv("../simulations/simulations3.csv")
  #d <- read_csv("simulations3.csv", col_types = cols())
  d <- d %>%
    mutate(Day = floor(Time / (60*60*24)) - 100,
           Hour = floor((d$Time / 3600) %% 24),
           RuminationFrequency=as_factor(RuminationFrequency))
  names(d)[18] <- "W"
  
  d <- d %>% 
    dplyr::select(Run, Day, PTEV, Gamma, PTES, W, 
                  NumAttributes, RuminationFrequency, 
                  MemoryEntropy, ChunkSimilarity, Traumatic, ChunkV) %>%
    filter(Day > -20)
  
  a <- d  %>%
    group_by(Run, Day, PTEV, PTES, Gamma, 
             NumAttributes, RuminationFrequency, W) %>%
    summarise(MemoryEntropy = mean(MemoryEntropy),
              ChunkSimilarity = mean(ChunkSimilarity),
              PTraumatic = mean(Traumatic),
              CTraumatic = sum(Traumatic),
              ChunkV = mean(ChunkV))
  
  write_csv(x=a, file = "../simulations/simulations3_aggregated2.csv")
}

a <- read_csv(gzfile("../simulations/simulations3_aggregated2.csv.gz"),
              col_types = cols())

Then, we are going to rename the simulation variables to make them consistent with the names used in published papers:

a <- a %>%
  rename(I = PTEV) %>%
  rename(gamma = Gamma) %>%
  rename(A = NumAttributes) %>%
  rename(C = PTES) %>%
  rename(R = RuminationFrequency)

Recovery Trajectories

To identify a recovery trajectory, we need to first define a classification function. Here, we are going to use the following algorithm:

First the model calculates three average values:

The mean probability P(baseline) of retrieving intrusive memories in the 10 days preceding the PTE, or baseline period. By definition, this is always zero.
The mean probability P(acute) of retrieving intrusive memories in the 10 days following the PTE, or during the acute period.
The mean probability P(chronic) of retrieving intrusive memories in the last ten days of the second month after the PTE (i.e, days 51-60), or the chronic period.

Then, the model calculates two statistical tests:

A t-test between the baseline and acute period, or acute test; and
A t-test between the acute-period and the chronic period, or chronic test.

And finally we apply the following decision tree:

If the acute test is significant at p < .05 (Pacute > Pbaseline), then

1.1. If the chronic test is also significant at p < .05, and Pchronic > Pacute, we classify the trajectory as Delayed.

1.2. If the chronic test is significant at p < .05 but P(chronic) < P(acute), then we classify the trajectory as Recovery;

1.3. If the chronic test is non-significant, then P(chronic) ~ P(acute) , and we classify the trajectory as Chronic.
If, instead, the acute test is not significant and P(acute) ~ P(baseline), then:

2.1. If the chronic test is significant at p < .05, then P(chronic) > P(acute), and we classify the trajectory as Delayed.

2.2 If the chronic test is also not significant, then P(chronic) ~ P(baseline) ~ P(acute), and we classify the trajectory as Resilient.

# Classifies people based on trajectory:
# Chronic = 1
# Recovery =2
# Delay    =3
# resilient= 4
#
classify <- function(days) {
  g1 <- days[10:19]
  g2 <- days[21:30]
  g3 <- days[70:79]
  t12 <- t.test(g1, g2)
  t23 <- t.test(g2, g3)
  category <- 0
  if (!is.na(t12$p.value) & t12$p.value < 0.05) {
    # t1 < t2
    if (!is.na(t23$p.value) & t23$p.value < 0.05) {
      # t2 <> t3

      if (!is.na(t23$statistic) & t23$statistic > 0) {
        # t1 < t2, t2 > 3
        category <- "Recovery" #2  # Recovery
      } else {
        # t1 < t2, t2 < t3
        category <- "Delayed" # 3 # Delayed
      }
    } else {
      # t1 < t2, t2 = t3
      category <- "Chronic" # 1 # Chronic
    }
  } else {
    # t1 == t2
    if (!is.na(t23$p.value) & t23$p.value < 0.05) {
      # t1 == t2, t2 < t3
      category <- "Delayed" #3 # Delayed
    } else {
      # t1 == t2, t2 = t3
      category <- "Resilient" # 4 # Resilient
    }
  }
  category
}

patterns <- a %>%
  group_by(Run, I, C, A, R, W, gamma) %>%
  dplyr::summarize(Trajectory = classify(PTraumatic))

## `summarise()` has grouped output by 'Run', 'I', 'C', 'A', 'R', 'W'. You can override using the `.groups` argument.

patterns <- patterns %>% ungroup

Finally, let’s assign each data point in the main averages with the corresponding trajectory, and re-order the trajectories in order of severity:

a <- inner_join(a, patterns)

## Joining, by = c("Run", "I", "C", "gamma", "A", "R", "W")

a$Trajectory <- factor(a$Trajectory,
                      levels = c("Recovery", "Delayed", "Resilient",  "Chronic"))

Chronic Probability of Traumatic Intrusions

The other “behavioral” measure is the incidence of traumatic memory intrusions. As an aggregate measure, we will consider only the averages in the “Chronic” period, which we will be taking as indicative of long-term consequences of the PTSD.

traumatic <- a %>%
  filter(Day > 1) %>%
  group_by(Run, I, C, R, W, gamma) %>%
  dplyr::summarize(PTraumatic = mean(PTraumatic),
            CTraumatic = sum(CTraumatic))

## `summarise()` has grouped output by 'Run', 'I', 'C', 'R', 'W'. You can override using the `.groups` argument.

At this point, we can visualize the mean aggregated timecourse of intrusions for each trajectory.

K = rev(brewer.pal(5, "Greys"))

TOTAL <- 28800  # Param combos

traj_total <- a %>%
  filter(Day == 30) %>%
  group_by(Trajectory) %>%
  summarise(N = length(CTraumatic),
            Prop = length(CTraumatic)/TOTAL,
            MeanC = mean(CTraumatic))

ggplot(filter(a, I != 1), 
       aes(x=Day, y=CTraumatic, col=Trajectory, fill=Trajectory)) +
  ggtitle("Memory Intrusions by Trajectory") +
  #geom_vline(data=props, aes(xintercept=-0.5)) +
  annotate("segment", x=-0, xend=-0,
           y=-Inf, yend=Inf, col="black", lty=1, size=1) +
  annotate("rect", xmin=50, xmax=60,
           ymin=-Inf, ymax=Inf, fill=K[1], alpha=0.3)+
  annotate("rect", xmin=1, xmax=10,
           ymin=-Inf, ymax=Inf, fill=K[3], alpha=0.3)+
  annotate("rect", xmin=-10, xmax=-1,
           ymin=-Inf, ymax=Inf, fill=K[4], alpha=0.3)+
  
  annotate("text", x= -5.5, y=15, label="Base")+
  annotate("text", x= 5.5, y=15, label="Acute")+
  annotate("text", x= 55, y=15, label="Chronic")+
  
  geom_text(data=traj_total, aes(x= 30, y = MeanC, 
                                 label=percent(Prop, accuracy = 0.1)),
            vjust=-.5, show.legend=F) +
  
  stat_summary(fun.data = mean_se, geom="line") +
  stat_summary(fun.data = mean_cl_normal, 
               geom="ribbon", alpha = 0.25, col=NA) +
  coord_cartesian(ylim=c(0, 15)) +
  ylab("Total Number of Memory Intrusions") +
  scale_color_d3() +
  scale_fill_d3() +
  theme_pander()

ggsave("figures/figure6.png", dpi = 300)

## Saving 5 x 3.5 in image

The results indicate that our classification was correct; the trajectories behave as planned.

Finding a baseline to compare the model to existing literature

Since we have no idea how typical the model parameters are of the existing patients and their cases, we need to provide some sort of baselines. To do so, we will use the proportion of trajectories that were given by Galatzer-Levy and Bonanno in their meta-analysis. They are:

Resilient, 0.657
Recovery, 0.208
Chronic, 0.106
Delayed, 0.089

Note that they do not normalize to one, due to the lack of some trajectories in some studies. To use them, we will need to make sure to normalized them first, which can be done at time of testing.

T_observed <- c(0.657, 0.208, 0.106, 0.089)
names(T_observed) <- c("Resilient", "Recovery", "Chronic", "Delayed")

We can now plot the percentages and the patterns by trajectory:

props <- patterns %>%
  group_by(Trajectory, 
           I, 
           gamma, 
           C, 
           W,
           R,
           #A,
           ) %>%
  dplyr::summarize(N = length(Trajectory)) %>%
  dplyr::mutate(Prop = N/50)

## `summarise()` has grouped output by 'Trajectory', 'I', 'gamma', 'C', 'W'. You can override using the `.groups` argument.

props %>%
  pivot_wider(id_cols=c(I, 
                        gamma, 
                        C,
                        W,
                        R,
                        #A,
                        ), 
              names_from = Trajectory,
              values_from = N,
              values_fill=0) -> test

mychi <- function(vals) {
  chisq.test(vals, 
             p=T_observed, 
             rescale.p=T)$statistic
}

mycorr <- function(vals) {
  cor(vals, T_observed)
}

myrmse <- function(vals) {
  v <- 50 * (T_observed/sum(T_observed))
  sqrt(mean((vals - v)^2))
}


lprops <- test %>%
  pivot_longer(cols=c("Resilient", "Recovery", "Chronic", "Delayed"),
               names_to = "Trajectory",
               values_to = "N")

lprops$Trajectory <- factor(lprops$Trajectory,
                      levels = c("Recovery", "Delayed", "Resilient",  "Chronic"))

# Suppresses X-squared approximation warnings
old.warn <- options()$warn
options(warn = -1)

atest <- lprops %>%
  group_by(I, 
           gamma, 
           C, 
           W,
           R,
           #A,
           ) %>%
  summarise(Chi = mychi(N), r=mycorr(N), RMSE=myrmse(N))

## `summarise()` has grouped output by 'I', 'gamma', 'C', 'W'. You can override using the `.groups` argument.

# resets warnings
options(warn = old.warn)

test <- inner_join(test, atest)

## Joining, by = c("I", "gamma", "C", "W", "R")

Now, we can identfy the baseline as the condition with the minimum \(\chi^2\)-squared test (or the minimum RMSE, they are the same)

This shows how correlation coefficient and \(\chi^2\) are related, while the RMSE does not really capture either of them:

ggplot(test, aes(x=Chi, y=r, col=RMSE)) + 
  geom_point(size=4, alpha=0.5) + 
  scale_color_viridis(option="plasma")+
  ylab(expression(italic(r))) +
  xlab(expression(chi^2)) +
  ggtitle("Relationship between Error Measures in Trajectories") +
  theme_pander()

Now, let’s identify a “baseline” condition—the one that most resembles the distribution of trajectories identified by Galatzer-Levy.

baseline <- test %>%
  filter(Chi == min(test$Chi))

baseline %>%
  kable(digits=3) %>%
  kable_styling(bootstrap_options = c("hover", "striped"))

I	gamma	C	W	R	Chronic	Delayed	Recovery	Resilient	Chi	r	RMSE
40	0.9	0.25	8	0	15	3	15	67	7.464	0.985	18.875

And here is a visual representation of the four trajectories at that baseline (smoothed, because there is too much day-to-day variance in this small sample).

baseline_a <- a %>%
  filter(I == baseline$I,
         W == baseline$W,
         C == baseline$C,
         R == baseline$R,
         gamma == baseline$gamma,
         )

base_total <- baseline_a %>%
  filter(Day == 15) %>%
  group_by(Trajectory) %>%
  summarise(N = length(CTraumatic),
            Prop = length(CTraumatic)/100,
            MeanC = mean(CTraumatic))

         
ggplot(baseline_a, 
       aes(x=Day, y=CTraumatic, col=Trajectory, fill=Trajectory)) +
  geom_smooth(method="loess", span=0.2) +
  ggtitle("Trajectories for Best-Fitting Parameters") +
  scale_color_d3() +
  geom_text(data=base_total, aes(x= c(25, 10, 20, 10), y = MeanC, 
                                 label=percent(Prop, accuracy = 0.1)),
            vjust=-.5, show.legend=F) +
  ylab("Total Number of Memory Intrusions") +
  geom_vline(data=props, aes(xintercept=-0.25)) +
  scale_fill_d3() +
  theme_pander()

## `geom_smooth()` using formula 'y ~ x'

ggsave("figures/figure7.png", dpi=300)

## Saving 5 x 3.5 in image
## `geom_smooth()` using formula 'y ~ x'

The Chonic and Delayed trajectories, being the least common, exhibit significant ups and downs, hence the need for Loess smoothing. Nonetheless, the trajectories remain visible and clearly identifiable.

But how good is the baseline condition, in terms of being closed to Galatzer-Levy’s estimated pooled values? We can just examine the significance of the \(\chi^2\) test.

baseline_vals <- baseline[c("Resilient", "Recovery", "Chronic", "Delayed")]
chi <- chisq.test(baseline_vals, p=T_observed, rescale.p = T)

The result is not too bad—A bit dissimilar, but not too far either, and not statistically significant (\(\chi^2\) = 7.463523, p = 0.058503)).

The main difference is that the delayed onset trajectory has a higher proportion than in the review. However, Galatzer-Levy and Bonanno themselves note that the resilient (M= 0.66) delayed onset trajectories (M = 0.099) has a higher prevalence in studies of single-point event. If these corrected prevalences are taken into account, the chi squared reduces further.

Factors affecting the distribution of trajectories

Although theoretically possible, a full five-way analysis of model’s factors and interactions would be statistically uninformative. Instead, we will examine the model’s outcomes and behaviors with the respect to the known main effects reported in the literature.

We can now examine, one by one, how the different factor in the model affect the distribution of trajectories.

To examine these factors, we need to use logistic regression. And, to do so, we need to create a number of dummy variables specifying a binary status for each trajectory in patterns.

patterns <- patterns %>%
  mutate(Resilient = if_else(Trajectory == "Resilient", 1, 0),
         Chronic = if_else(Trajectory == "Chronic", 1, 0),
         Recovery = if_else(Trajectory == "Recovery", 1, 0),
         Delayed = if_else(Trajectory == "Delayed", 1, 0))

External factors: PTE Intensity and Context Similarity

There are two main external factors in the model: The intensity of the PTE (I) and the contextual similarity (C). They will be both examined here.

Effects on Timecourses

First, let’s load additional values of I and C in addition to the ones in the baseline.

IC_a <- a %>%
  filter(#I == baseline$I,
    I != 1,
         W == baseline$W,
         #C == baseline$C,
         R == baseline$R,
         gamma == baseline$gamma,
         #A == baseline$A
         )

IC_patterns <- patterns %>%
    filter(#I == baseline$I,
    I != 1,
         W == baseline$W,
         #C == baseline$C,
         R == baseline$R,
         gamma == baseline$gamma,
         #A == baseline$A
         )

IC_trajectories <- lprops %>%
    filter(#I == baseline$I,
    I != 1,
         W == baseline$W,
         #C == baseline$C,
         R == baseline$R,
         gamma == baseline$gamma,
         #A == baseline$A
         )


IC_traumatic <- traumatic %>%
      filter(#I == baseline$I,
         I != 1,
         W == as.numeric(baseline$W),
         #C == baseline$C,
         R == as.numeric(baseline$R),
         gamma == as.numeric(baseline$gamma),
         #A == baseline$A
         )

First, let’s look at the time course of intrusions:

ggplot(IC_a, 
       aes(x=Day, y=CTraumatic, col=as.factor(I), fill=as.factor(I))) +
  #geom_smooth(method="loess", span=0.2) +
  stat_summary(fun.data = mean_se, geom="line") +
  stat_summary(fun.data = mean_cl_normal, 
               geom="ribbon", alpha = 0.25, col=NA) +
  
  facet_wrap(~C, labeller = label_bquote(italic(C) == .(C))) + #label=label_both) +
  labs(col=expression(italic(I)[PTE]), 
       fill=expression(italic(I)[PTE]), 
       facet="Z") +
  #coord_cartesian(ylim=c(0, 0.5)) +
  ggtitle("Effects of PTE Intensity") +
  ylab("Number of Traumatic Memory Intrusions") +
  
  scale_color_viridis(discrete=T, option="plasma", end=0.8) +
  scale_fill_viridis(option = "plasma", discrete = T, end=0.8) +
  geom_vline(data=props, aes(xintercept=-0.15)) +
  theme_pander()

The effects of I and C on the probabilities of recall during the Chronic period can be understood with a linear model analysis

IC_mod <- glm(CTraumatic ~ I * C, family="poisson", IC_traumatic)
summary(IC_mod) %>%
  xtable %>%
  kable %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

	Estimate	Std. Error	z value
(Intercept)	-0.7993158	0.0439928	-18.16923
I	0.0842682	0.0008007	105.23676
C	7.8729269	0.0636742	123.64387
I:C	-0.0739858	0.0011709	-63.18745

Effects on Trajectories

And here are the changes in trajectories.

ggplot(data=IC_trajectories, aes(x="", y=N/100, fill=Trajectory)) +
  geom_bar(stat = "identity", col="grey", width=1) +
  #facet_grid(C~I, labeller = label_both) +
  facet_grid(C~I, labeller = label_bquote(
    rows = italic(C) == .(C),
    cols = italic(I)[PTE] == .(I))) +
  coord_polar("y", start=0) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 2L)) +
  scale_fill_d3() +
  scale_color_d3() +
  ylab("PTE Intensity") +
  xlab("Contextual Similarity") +
  ggtitle("Effect of External Factors on Trajectories") +
  geom_text_repel(aes(label=percent(N/100, .1)), col="white",
            position=position_stack(vjust=0.5 ), size=3)+
  theme_pander() +
  theme(axis.ticks = element_blank(), 
        axis.text.y = element_blank(),
        axis.text.x = element_blank(),
        legend.position = "bottom")

Are the change in trajectories meaningful? To do so, we first examine the trajectories using a MANOVA:

man <- manova(cbind(Resilient, Recovery, Delayed) ~ I * C, IC_patterns)
tidy(man) %>%
  xtable() %>%
  kable() %>%
  kable_styling(bootstrap_options = c("hover", "striped"))

term	df	pillai	statistic	num.df	den.df	p.value
I	1	0.2505389	133.0483	3	1194	0
C	1	0.5364246	460.5442	3	1194	0
I:C	1	0.0589542	24.9337	3	1194	0
Residuals	1196	NA	NA	NA	NA	NA

Then, we examine each category independently, we will perform a logistic regression on dummy binary variables that categorize each run of the model as belonging to that category, or not.

In the case of the Resilient and Recovery trajectories, both intensity I, Context C are significant, but their interaction is not. This implies that the effects of I and C are additive.

Furthermore, for the Resilient trajectories, the coefficients are negative, implying that the proportion of Recovery and Resilient trajectories decreases as the intensity and the context similarity increase.

Resilient Trajectory

mylogit <- glm(Resilient ~ (I*C), family="binomial", data=IC_patterns)
summary(mylogit) %>%
  xtable() %>%
  kable() %>%
  kable_styling(bootstrap_options = c("hover", "striped"))

	Estimate	Std. Error	z value	Pr(>\|z\|)
(Intercept)	8.0726911	0.7396808	10.9137497	0.0000000
I	-0.1042898	0.0141387	-7.3761937	0.0000000
C	-10.3446128	1.4043857	-7.3659343	0.0000000
I:C	-0.0144672	0.0322948	-0.4479746	0.6541716

Recovery Trajectory

For Recovery trajectories, however, the coefficients are significantly positive, implying that the number of cases increase with the severity of intensity and similarity. This is likely due to a greater number of resilient models developing intrusive memories following PTE under greater contextual similarity.

mylogit <- glm(Recovery ~ (I*C), family="binomial", data=IC_patterns)
summary(mylogit) %>%
  xtable() %>%
  kable() %>%
  kable_styling(bootstrap_options = c("hover", "striped"))

	Estimate	Std. Error	z value	Pr(>\|z\|)
(Intercept)	-6.5266334	0.7852632	-8.3113963	0.0000000
I	0.0757766	0.0147198	5.1479313	0.0000003
C	5.0969156	1.2604091	4.0438580	0.0000526
I:C	-0.0085358	0.0247694	-0.3446121	0.7303860

Chronic Trajector

In the case of Chronic and Delayed trajectories, however, the pattern is different. In these cases, both the main factors and their interaction are significant. Furthermore, in both cases, the main factors have positive coefficients but their interaction is negative. This suggests that the effects of one factor decrease the size of the other factor increases.

mylogit <- glm(Chronic ~ (I*C), family="binomial", data=IC_patterns)
summary(mylogit) %>%
  xtable() %>%
  kable() %>%
  kable_styling(bootstrap_options = c("hover", "striped"))

	Estimate	Std. Error	z value
(Intercept)	-6.2050051	0.6318950	-9.819677
I	0.0760648	0.0121663	6.252087
C	9.4454215	1.0400816	9.081424
I:C	-0.1332185	0.0207257	-6.427706

Delayed Trajectory

mylogit <- glm(Delayed ~ (I*C), family="binomial", data=IC_patterns)
summary(mylogit) %>%
  xtable() %>%
  kable() %>%
  kable_styling(bootstrap_options = c("hover", "striped"))

	Estimate	Std. Error	z value	Pr(>\|z\|)
(Intercept)	-9.2380905	1.7014060	-5.429680	0.0000001
I	0.1004351	0.0311299	3.226324	0.0012539
C	11.4814812	2.4891026	4.612699	0.0000040
I:C	-0.2098212	0.0489880	-4.283111	0.0000184

Internal Factors: Working Memory and Recollection Vividness

The internal factors are two, differences in WMC and differences in vividness.

This can be related to the difference in children.

Wg_a <- a %>%
  filter(I == baseline$I,
         #W == baseline$W,
         C == baseline$C,
         R == baseline$R,
         #gamma == baseline$gamma,
         #A == baseline$A
         )

Wg_test <- lprops %>%
    filter(I == baseline$I,
         #W == baseline$W,
         C == baseline$C,
         R == baseline$R,
         #gamma == baseline$gamma,
         #A == baseline$A
         )

Wg_patterns <- patterns %>%
    filter(I == baseline$I,
         #I != 1,
         #W == as.numeric(baseline$W),
         C == baseline$C,
         R == as.numeric(baseline$R),
         #gamma == as.numeric(baseline$gamma),
         #A == baseline$A
         )

Wg_traumatic <- traumatic %>%
    filter(I == baseline$I,
         #I != 1,
         #W == as.numeric(baseline$W),
         C == baseline$C,
         R == as.numeric(baseline$R),
         #gamma == as.numeric(baseline$gamma),
         #A == baseline$A
         )

Effects on the Timecourse of Intrusive Memories

Here we visualize the change in timecourses.

Wg_a$gamma <- as_factor(Wg_a$gamma)
ggplot(Wg_a, 
       aes(x=Day, y=CTraumatic, col=gamma, fill=gamma)) +
#  geom_point(alpha=0.1) +
  stat_summary(fun.data = mean_se, geom="line") +
  stat_summary(fun.data = mean_cl_normal, 
               geom="ribbon", alpha = 0.25, col=NA) +
  
  #geom_smooth(method="loess", span=0.2) +
  facet_wrap(~W, ncol=3,
             labeller  = label_bquote(
               rows = italic(W) == .(W))) +
  geom_vline(data=props, aes(xintercept=-0.15)) +
  ylab("Number of\nTraumatic Memory Intrusions") +
  #coord_cartesian(ylim=c(0, 0.75)) +
  ggtitle("Effects of Idiographic Parameters") +
  scale_color_viridis(expression(paste(gamma, " =")), discrete=T, 
                      option="plasma", begin = 0, end=0.9) +
  scale_fill_viridis(expression(paste(gamma, " =")), option = "plasma", 
                     discrete = T, begin = 0, end=0.9) +
  theme_pander() +
  theme(legend.position = "bottom")

It is striking to notice that, within baseline levels of I and C, the effects of working memory are incredibly large. This is likely due to the fact that the parameters were not properly calibrated.

And here is the corresponding analysis. Like before, the main factors W and \(\gamma\) are significant. However, this time, their interaction is not, suggesting that their effects are additive.

Wg_mod <- glm(CTraumatic ~ W * gamma, family="poisson", Wg_traumatic)
summary(Wg_mod) %>%
  xtable %>%
  kable %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

	Estimate	Std. Error	z value	Pr(>\|z\|)
(Intercept)	-4.6088943	0.2487333	-18.5294648	0.0000000
W	-0.5902732	0.0547969	-10.7720158	0.0000000
gamma	15.0052321	0.2687441	55.8346535	0.0000000
W:gamma	-0.0295063	0.0592083	-0.4983467	0.6182397

Effects on Trajectories

And here are the changes in trajectories.

ggplot(data=Wg_test, aes(x="", y=N/100, fill=Trajectory)) +
  geom_bar(stat = "identity", col="white", width=1, alpha=1) +
  facet_grid(gamma~W, 
             labeller  = label_bquote(
               rows = gamma == .(gamma),
               cols = italic(W) == .(W))) +
  coord_polar("y", start=0) +
  #scale_fill_d3() +
  scale_y_continuous(labels = scales::percent_format(accuracy = 2L)) +
  scale_fill_d3() +
  xlab("Recollection Vividness") +
  ylab("Cognitive Control") +
  ggtitle("Effect of Idiographic Factors on Trajectories") +
  geom_text_repel(aes(label=percent(N/100, .1)),
            position=position_stack(vjust=0.5), size=3, col="white") +
  #geom_text_repel() +
  theme_pander() +
  theme(axis.ticks = element_blank(), 
        axis.text.y = element_blank(),
        axis.text.x = element_blank(),
        legend.position = "bottom")

As before, we can start with a MANOVA:

man <- manova(cbind(Resilient, Recovery, Delayed) ~ W * gamma, Wg_patterns)
tidy(man) %>%
  xtable() %>%
  kable() %>%
  kable_styling(bootstrap_options = c("hover", "striped"))

term	df	pillai	statistic	num.df	den.df	p.value
W	1	0.6151552	476.33812	3	894	0
gamma	1	0.0818391	26.56184	3	894	0
W:gamma	1	0.0478437	14.97384	3	894	0
Residuals	896	NA	NA	NA	NA	NA

The MANOVA shows that both factors and their interactions are significant. We can now carry out a trajectory by trajectory analysis.

Because of the large varince in the proportion of trajectories (due to the massive effects of W), none of the interaction terms are significant. Thus, we analyzed all four trajectories with a simpler generalized linear model of the form \(Y = \beta_0 + \beta_W W + \beta_{\gamma} \gamma + \epsilon\), using a logistic link function.

Resilient Trajectory

mylogit <- glm(Resilient ~ W + gamma, 
                  family = binomial("logit"),
                  data = Wg_patterns)
summary(mylogit)  %>%
  xtable() %>%
  kable() %>%
  kable_styling(bootstrap_options = c("hover", "striped"))

	Estimate	Std. Error	z value	Pr(>\|z\|)
(Intercept)	4.751373	1.771695	2.681823	0.0073222
W	1.034876	0.071035	14.568530	0.0000000
gamma	-13.553966	2.178413	-6.221945	0.0000000

Surprisingly, there is no significant effect of either parameter on the Resilient trajectory. This is likely due to the massive amount of variance, with the Resilient Trajectory going from being almost absent when W = 4 to being the only possible trajectory when W = 12.

Chronic Trajectory

We can then examine the Chronic trajectory:

mylogit <- glm(Chronic ~ (W + gamma), family=binomial("logit"), 
               data=Wg_patterns)
summary(mylogit) %>%
  xtable() %>%
  kable() %>%
  kable_styling(bootstrap_options = c("hover", "striped"))

	Estimate	Std. Error	z value	Pr(>\|z\|)
(Intercept)	-6.5273256	1.5240563	-4.282864	1.85e-05
W	-0.4288011	0.0402964	-10.641173	0.00e+00
gamma	8.8505973	1.7098097	5.176364	2.00e-07

Recovery Trajectory

Then, we on the Recovery trajectory.

mylogit <- glm(Recovery ~ (W + gamma), family=binomial("logit"), 
               data=filter(Wg_patterns, W<13))
summary(mylogit) %>%
  xtable() %>%
  kable() %>%
  kable_styling(bootstrap_options = c("hover", "striped"))

	Estimate	Std. Error	z value	Pr(>\|z\|)
(Intercept)	7.0495506	1.5366820	4.587514	0.0000045
W	-0.6371039	0.0585618	-10.879166	0.0000000
gamma	-5.3828643	1.6813203	-3.201570	0.0013668

Delayed Trajectory

The Delayed trajectory:

mylogit <- glm(Delayed ~ (W + gamma), 
               family=binomial("logit"), 
               data=filter(Wg_patterns, W < 20) )
summary(mylogit) %>%
  xtable() %>%
  kable() %>%
  kable_styling(bootstrap_options = c("hover", "striped"))

	Estimate	Std. Error	z value	Pr(>\|z\|)
(Intercept)	-13.3386513	2.9206263	-4.567052	4.9e-06
W	-0.3133589	0.0566509	-5.531406	0.0e+00
gamma	14.1075665	3.1741266	4.444551	8.8e-06

Comorbidty: Effects of Rumination

R_a <- a %>%
  filter(I == baseline$I,
         W == baseline$W,
         C == baseline$C,
         #R == baseline$R,
         gamma == baseline$gamma,
         #A == baseline$A
         )

R_test <- lprops %>%
    filter(I == baseline$I,
         W == baseline$W,
         C == baseline$C,
         #R == baseline$R,
         gamma == baseline$gamma,
         #A == baseline$A
         )

R_patterns <- patterns %>%
    filter(I == baseline$I,
         #I != 1,
         W == as.numeric(baseline$W),
         C == baseline$C,
         #R == as.numeric(baseline$R),
         gamma == as.numeric(baseline$gamma),
         #A == baseline$A
         )

R_traumatic <- traumatic %>%
    filter(I == baseline$I,
         #I != 1,
         W == as.numeric(baseline$W),
         C == baseline$C,
         #R == as.numeric(baseline$R),
         gamma == as.numeric(baseline$gamma),
         #A == baseline$A
         )

Effects on the Timecourse of Intrusive Memories

Here we visualize the change in timecourses.

ggplot(R_a, 
       aes(x=Day, y=CTraumatic, col=as.factor(R), fill=as.factor(R))) +
  stat_summary(fun.data = mean_se, geom="line") +
  stat_summary(fun.data = mean_cl_normal, 
               geom="ribbon", alpha = 0.25, col=NA) +
  geom_vline(data=props, aes(xintercept=-0.15)) +

  #geom_smooth(method="loess", span=0.2) +
  #facet_wrap(~R) +
  #coord_cartesian(ylim=c(0, 0.75)) +
  ylab("Number of\nTraumatic Memory Intrusions") +
  #coord_cartesian(ylim=c(0, 0.75)) +
  ggtitle("Effects of Rumination") +
  scale_color_viridis(expression(paste(italic(R), " =")), discrete=T, option="plasma", end=0.8) +
  scale_fill_viridis(expression(paste(italic(R), " =")), option = "plasma", discrete = T, end=0.8) +
  theme_pander() +
  theme(legend.position = "bottom")

The effects of Rumination are huge. This is likely due to the fact that the parameters were not properly calibrated.

And here is the corresponding analysis. Like before, the main factors W and \(\gamma\) are significant. However, this time, their interaction is not, suggesting that their effects are additive.

R_mod <- glm(CTraumatic ~ R, family="poisson", R_traumatic)
summary(R_mod) %>%
  xtable %>%
  kable %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

	Estimate	Std. Error	z value	Pr(>\|z\|)
(Intercept)	4.0167433	0.0189797	211.6333	0
R	0.1350059	0.0009804	137.7115	0

Effects on Trajectories

And here are the changes in trajectories.

ggplot(data=R_test, aes(x="", y=N/100, fill=Trajectory)) +
  geom_bar(stat = "identity", col="white", width=1, alpha=1) +
  facet_wrap(~ R, labeller = label_bquote(
               rows = italic(R) == .(R))) +
  coord_polar("y", start=0) +
  #scale_fill_d3() +
  scale_y_continuous(labels = scales::percent_format(accuracy = 2L)) +
  scale_fill_d3() +
  ylab("Rumination") +
  xlab("") +
  ggtitle("Effect of Rumination on Trajectories") +
  geom_text_repel(aes(label=percent(N/100, .1)),
            position=position_stack(vjust=0.5), size=3, col="white") +
  #geom_text_repel() +
  theme_pander() +
  theme(axis.ticks = element_blank(), 
        axis.text.y = element_blank(),
        axis.text.x = element_blank())

As before, we can start with a MANOVA:

man <- manova(cbind(Resilient, Recovery, Delayed) ~ R, R_patterns)
tidy(man) %>%
  xtable() %>%
  kable() %>%
  kable_styling(bootstrap_options = c("hover", "striped"))

term	df	pillai	statistic	num.df	den.df	p.value
R	1	0.4574771	55.09169	3	196	0
Residuals	198	NA	NA	NA	NA	NA

The MANOVA shows that rumination is significant. We can now carry out a trajectory by trajectory analysis.

Resilient Trajectory

First, we can examine the effects on the Resilient trajectory.

mylogit <- glm(Resilient ~ R, 
               family = "binomial",
               data = R_patterns)
summary(mylogit) %>%
  xtable() %>%
  kable() %>%
  kable_styling(bootstrap_options = c("hover", "striped"))

	Estimate	Std. Error	z value	Pr(>\|z\|)
(Intercept)	0.7081851	0.2126697	3.329976	0.0008685
R	-0.2092142	0.0311797	-6.709947	0.0000000

Chronic Trajectory

We can then examine the Chronic trajectory:

mylogit <- glm(Chronic ~ R, family="binomial", data=R_patterns)
summary(mylogit) %>%
  xtable() %>%
  kable() %>%
  kable_styling(bootstrap_options = c("hover", "striped"))

	Estimate	Std. Error	z value	Pr(>\|z\|)
(Intercept)	-1.7346011	0.2800560	-6.193765	0.0000000
R	0.0622526	0.0173836	3.581108	0.0003421

Recovery Trajctory

The, we analyze the Recovery trajectory.

mylogit <- glm(Recovery ~ R, family="binomial", data = R_patterns)
summary(mylogit) %>%
  xtable() %>%
  kable() %>%
  kable_styling(bootstrap_options = c("hover", "striped"))

	Estimate	Std. Error	z value	Pr(>\|z\|)
(Intercept)	-1.7346010	0.2800209	-6.194541	0
R	0.0967636	0.0172348	5.614430	0

And, finally, The Delayed trajectory:

mylogit <- glm(Delayed ~ R, 
               family="binomial", 
               data=R_patterns)
summary(mylogit) %>%
  xtable() %>%
  kable() %>%
  kable_styling(bootstrap_options = c("hover", "striped"))

	Estimate	Std. Error	z value	Pr(>\|z\|)
(Intercept)	-3.4760987	0.5862093	-5.9297915	0.0000000
R	0.0149022	0.0388606	0.3834797	0.7013641

Neurobiology of PTSD

Finally, we will examine how the model captures some known effects of the neurobiology of PTSD and, in particular, the size of the hippocampus.

Hippocampus Size

First, we define a baseline condition: No PTE (I = 1) and \(\gamma\) = 0.8.

H1 <- a %>%
  filter(I == 1 & gamma == 0.9 & C==0.25 & R == 0 & Day > 49) %>%
  dplyr::select(MemoryEntropy) %>%
  summarise(H=mean(MemoryEntropy)) %>%
  as.double()

Then, we define a function that calculates memory entropy on the last N days of a timeseries.

h_entropy <- function(timeseries, N = 10) {
  L <- length(timeseries)
  mean(timeseries[L-N:N])
}

With this in place, we can now calculate the predicted decrease in hippocampus size.

hsize <- a %>%
  filter(Day > 49) %>%
  group_by(Run, I, C, R, W, gamma, Trajectory) %>%
  dplyr::summarize(H = h_entropy(MemoryEntropy) - H1,
            HippocampusDecrease = 100*(mean(MemoryEntropy)/H1 - 1),
            MeanTraumatic = mean(PTraumatic),
            NumTraumatic = sum(CTraumatic),
            MeanSimilarity = mean(ChunkSimilarity),
            Trajectory = )

## `summarise()` has grouped output by 'Run', 'I', 'C', 'R', 'W', 'gamma'. You can override using the `.groups` argument.

Visualizing Hippocampus Volume Decrease

For pretty labels, we define two new labeller functions:

ptev_val <- function(val) {
  paste("I =", val)
}

gamma_val <- function(val) {
  #expression(paste(gamma, "=", val))
  paste("gamma =", val)
}

hsize <- hsize %>%
  mutate(External = paste("I = ", I , ", C = ", C))

ggplot(filter(hsize, I != 1),
       aes(x=NumTraumatic, y=HippocampusDecrease, col = Trajectory)) +
  geom_point(alpha=0.1, size=0.1, pch=21) +
  #geom_density_2d() +
  facet_grid( W ~ gamma, labeller=label_both) +
  #scale_color_brewer(palette="Paired") +
  scale_color_viridis(option = "plasma", discrete = T, end = 0.8) +
  stat_summary(fun.data = mean_se, geom="point", size=1) +
  geom_smooth(method = "lm", formula = y ~ x, col="black", fullrange= F) +
  #geom_smooth(method = "lm", aes(fill=Condition), se=F) +
  ylab("Percentage Decrease in Hippocampus Volume") +
  xlab("Cumulative Number of Traumatic Memory Intrusions (Day 50-60)") +
  ggtitle("Symptom Severity And Hippocampal Volume Decrease") +
  theme_pander() +
#  theme(legend.position = "bottom") +
  theme(panel.background=element_rect(fill="NA", colour="black"))

By Trajectory

hsize_trajectory <- a %>%
  filter(Day > 49, I != 1) %>%
  group_by(Run, I, C, R, W, gamma, Trajectory) %>%
  dplyr::summarize(HippocampusDecrease = 100*(mean(MemoryEntropy)/H1 - 1),
                   PTraumatic = mean(PTraumatic),
                   CTraumatic = sum(CTraumatic),
                   MeanSimilarity = mean(ChunkSimilarity))

## `summarise()` has grouped output by 'Run', 'I', 'C', 'R', 'W', 'gamma'. You can override using the `.groups` argument.

And here is the corresponding plot:

ggplot(hsize_trajectory, aes(x=Trajectory, y = HippocampusDecrease, 
                             col=Trajectory, fill=Trajectory)) +
  facet_wrap(~I, labeller = label_bquote(
    rows = italic(I)[PTE] == .(I))) +
  scale_fill_d3() +
  scale_color_d3() +
  stat_summary(geom="point", fun.data="mean_sdl") +
  stat_summary(geom="errorbar", fun.data="mean_sdl", width=.25, 
               fun.args = list(mult=1)) +
  #geom_boxplot(width=.25, fill=NA) +
  geom_violin(alpha=0.4, col="NA", trim = T, width  = 1.5) +
  theme_pander() +
  ylab("% Hippocampus Volume Decrease") +
  theme(axis.text.x = element_text(angle = 90)) +
  theme(panel.background=element_rect(fill="NA", colour="black"),
        legend.position = "bottom")

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

In general, largest decreases of hippocampus occur in the Chronic trajectory, and the smaller decreases in the resilient ones. But is a larger decrease inverse;y associated with the proportion of resilient trajectories?

hsize_traj_prop <- hsize_trajectory %>%
  mutate(IsResilient = ifelse(Trajectory=="Resilient", 1, 0)) %>%
  group_by(I, C, R, W, gamma) %>%
  summarise(HippocampusDecrease = mean(HippocampusDecrease),
          PropResilient = mean(IsResilient),
          LogVolume = log(100 + HippocampusDecrease))

## `summarise()` has grouped output by 'I', 'C', 'R', 'W'. You can override using the `.groups` argument.

Now, we can also replicate Briana’s original analysis of the correlation between PTSD and symptoms, but coloring the points on the basis of teir trajectory.

ggplot(filter(hsize_trajectory, I != 1),
       aes(x=CTraumatic, y=HippocampusDecrease, col = Trajectory)) +
  geom_point(alpha = .2, size=0.5) +
  #geom_density_2d() +
  scale_color_d3() + 
  #scale_color_brewer(palette="Dark2") +
  #scale_color_viridis(discrete=T, option="inferno", end=0.8) +
  #stat_summary(fun.data = mean_se, geom="point", size=1) +
  #geom_smooth(method = "lm", formula = y ~ x, aes(col=as.factor(I)), fill="white",
  #            fullrange= T) + # theme_pander() +
  geom_smooth(method = "lm", aes(fill=Trajectory), 
              se=T, line_style="dashed") +
  scale_fill_d3() + 
  ylab("Percentage Decrease in Hippocampus Volume") +
  xlab("Total Number of Memory Intrusions (Day 50-60)") +
  ggtitle("Correlation Between Symptom Severity\nAnd Hippocampal Volume Decrease") +
  theme_pander() +
  theme(legend.position = "bottom") #+

## Warning: Ignoring unknown parameters: line_style

## `geom_smooth()` using formula 'y ~ x'

Finally, we can now visualize the data on the DK brain atlas for clarity.

hsize_trajectory$Group = "PTSD"
hsize_trajectory$Group[hsize_trajectory$Trajectory == "Resilient"] <- "Control"

hpc_group <- hsize_trajectory %>%
  group_by(Group, Trajectory) %>%
  summarise(VolumeChange = mean(HippocampusDecrease),
            SD = sd(HippocampusDecrease)) %>%
  add_column(region = "parahippocampal")

## `summarise()` has grouped output by 'Group'. You can override using the `.groups` argument.

hpc_group$Trajectory <- as.character(hpc_group$Trajectory)

p1 <- hpc_group %>%
  group_by(Group, Trajectory) %>%
  ggplot(aes(fill = VolumeChange)) +
  geom_brain(atlas = dk, 
             col="white",
             show.legend = T) +
  coord_sf(xlim = c(750, 1050)) +
  facet_wrap(~ Group + Trajectory, 
             ncol=4,
             labeller = label_bquote(
               rows =  .(Group)~( .(Trajectory)))) +
  scale_fill_viridis("Mean %\nVolume Change", 
                     option = "plasma", 
                     na.value="lightgrey",
                     direction = 1,
                     limits=c(-10.0, 0)
                     ) +
  theme_brain2(text.family = "sans",
               text.colour="black") +
  theme(strip.background = element_blank(),
        strip.text = element_text(size=10, face="bold"),
        legend.position = "right",
        legend.title = element_text(size=10))

p2 <- ggplot(hsize_trajectory, 
             aes(x = "",
                 y = HippocampusDecrease)) +
  geom_violin(aes(fill=Group), alpha=0.5, 
              col="white", width=1.05, bw=2) +
  scale_fill_npg() +
  stat_summary(geom="point", fun.data="mean_sdl") +
  facet_wrap(~ Group + Trajectory, ncol=4) +
  
  stat_summary(geom="errorbar", fun.data="mean_sdl", width=.1, 
               fun.args = list(mult=1)) +
  ylab("% Volume Change") +
  xlab("") +
  theme_minimal() +
  theme(strip.text=element_blank())

p1 / p2

## merging atlas and data by 'region'

## Warning in strip_mat[panel_pos] <- unlist(unname(strips), recursive = FALSE)
## [[params$strip.position]]: number of items to replace is not a multiple of
## replacement length

ggsave("figures/figure15.png", dpi=300)

## Saving 8 x 4 in image

## merging atlas and data by 'region'

## Warning in strip_mat[panel_pos] <- unlist(unname(strips), recursive = FALSE)
## [[params$strip.position]]: number of items to replace is not a multiple of
## replacement length

Comparisons

Let’s look at some raw data

hsize_trajectory %>%
  group_by(Group) %>%
  summarise(VolumeChange = mean(HippocampusDecrease),
            SD = sd(HippocampusDecrease)) %>%
  xtable() %>%
  kable() %>%
  kable_styling(bootstrap_options=c("striped", "hover"))

Group	VolumeChange	SD
Control	-0.2626291	1.688191
PTSD	-9.0538188	7.771086

As expected, the the difference between simulated control groups (Traumatized controls and Traumatized PTSD) is significant to a two-sided, unequal variance t-test.

t.test(HippocampusDecrease ~ Group, 
       hsize_trajectory, 
       var.equal=F) %>%
  tidy() %>%
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

estimate	estimate1	estimate2	statistic	p.value	parameter	conf.low	conf.high	method	alternative
8.79119	-0.2626291	-9.053819	110.8244	0	12754.24	8.6357	8.946679	Welch Two Sample t-test	two.sided

Also, the the Control group shows a decrease of hippocampal volume that is small but significantly different than zero:

t.test(filter(hsize_trajectory, Group == "Control")$HippocampusDecrease,
       mean = 0) %>%
  tidy() %>%
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

estimate	statistic	p.value	parameter	conf.low	conf.high	method	alternative
-0.2626291	-10.53394	0	4584	-0.3115072	-0.2137509	One Sample t-test	two.sided

Finally, in the PTSD group, there is a significant difference by trajectories.

traumatized_group_hsize <- hsize_trajectory %>%
  filter(Group == "PTSD")

traumatized_group_hsize %>%
  aov(HippocampusDecrease ~ Trajectory, .)%>% 
  summary() %>%
  xtable() %>% 
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

	Df	Sum Sq	Mean Sq	F value	Pr(>F)
Trajectory	2	13121.93	6560.96562	110.886	0
Residuals	10646	629908.36	59.16855	NA	NA

A pairwise comparison indicates which trajectories are different:

pairwise.t.test(traumatized_group_hsize$HippocampusDecrease, 
                traumatized_group_hsize$Trajectory, 
                p.adjust.method = "bonferroni", paired=F, pool.sd = F) %>%
  tidy() %>%
  kable() %>%
  kable_styling(bootstrap_options = c("stiped", "hover"))

group1	group2	p.value
Delayed	Recovery	0
Chronic	Recovery	0
Chronic	Delayed	0

Functional Connectivity

Functional conectivity is measured through the similarity between each retrieved chunk and the current context, which is recorded in the ChunkSimilarity variable.

We can quickly visualize the effects of different trajectories on connectivity.

k <- a %>%
  filter(! is.na(ChunkSimilarity),
         C == 0,
         Day > 50)

ggplot(filter(k, I > 1), 
       aes(x=Trajectory, y = ChunkSimilarity, 
           col=Trajectory, fill=Trajectory)) +
  stat_summary(fun.data=mean_se, geom="bar", alpha=0.5) +
  stat_summary(fun.data=mean_cl_boot, geom="errorbar", width=0.25, size=1) +
  facet_wrap(~I, ncol=3, 
             labeller  = label_bquote(
               rows = italic(I)[PTE] == .(I))) +
  ylab("PFC - Hippocampus Connectivity") +
  scale_color_d3() +
  scale_fill_d3() +
  theme_pander() +
  theme(axis.text.x = element_blank(),
        legend.position = "bottom")

Visualizing Functional Connectivity

And now, we can project the mean differences in connectivity onto the DK brain atlas. The analysis is performed by simulated patient group, symptomatic PTSD (the set of Recovery, Chronic, and Delayed trajectories) vs. asymptomatic traumatized group (Resilient trajectory).

aconn <- a %>%
  filter(! is.na(ChunkSimilarity),
         C == 0,
         I > 1,
         Day > 50)

aconn$Trajectory <- as.character(aconn$Trajectory)

aconn$Group <- "PTSD"
aconn$Group[aconn$Trajectory == "Resilient"] <- "Control"

conn_group <- aconn %>%
  group_by(Group, Trajectory) %>%
  summarise(FC = mean(ChunkSimilarity),
            SD = sd(ChunkSimilarity))

## `summarise()` has grouped output by 'Group'. You can override using the `.groups` argument.

conn1 <- conn_group %>%
  add_column(region = "parahippocampal") 

conn2 <- conn_group %>%
  add_column(region = "medial orbitofrontal")

conn_group <- rbind(conn1, conn2)

p1 <- conn_group %>%
  group_by(Group, Trajectory) %>%
  ggplot(aes(fill = FC)) +
  geom_brain(atlas = dk, 
             col="white",
             show.legend = T) +
  coord_sf(xlim = c(750, 1050)) +
  facet_wrap(~ Group + Trajectory, 
             ncol=4,
             labeller = label_bquote(
               rows =  .(Group)~( .(Trajectory)))) +
  scale_fill_viridis("Mean Functional\nConnectivity", 
                     option = "plasma", 
                     na.value="lightgrey",
                     direction = 1,
                     limits=c(0.0, 0.5)
                     ) +
  theme_brain2(text.family = "sans",
               text.colour="black") +
  theme(strip.background = element_blank(),
        strip.text = element_text(size=10, face="bold"),
        legend.position = "right",
        legend.title = element_text(size=10))

p2 <- ggplot(aconn, 
             aes(x = "",
                 y = ChunkSimilarity)) +
  geom_violin(aes(fill=Group), alpha=0.5, 
              col="white", width=1.05,
              bw=0.075) +
  scale_fill_npg() +
  stat_summary(geom="point", fun.data="mean_sdl") +
  stat_summary(geom="errorbar", fun.data="mean_sdl", width=.1, 
               fun.args = list(mult=1)) +
  facet_wrap(~ Group + Trajectory, ncol=4) +
  ylab("Functional Connectivity") +
  xlab("") +
  theme_minimal() +
  theme(strip.text=element_blank())

# Multi-plot
p1 / p2

## merging atlas and data by 'region'

## Warning in strip_mat[panel_pos] <- unlist(unname(strips), recursive = FALSE)
## [[params$strip.position]]: number of items to replace is not a multiple of
## replacement length

ggsave("figures/figure16.png", dpi=300)

## Saving 8 x 4 in image

## merging atlas and data by 'region'

## Warning in strip_mat[panel_pos] <- unlist(unname(strips), recursive = FALSE)
## [[params$strip.position]]: number of items to replace is not a multiple of
## replacement length

Comparisons

Let’s look at the raw data and some descriptive statistics about predicted functional connectivity:

conn_group %>%
  filter(region == "parahippocampal") %>%
  select(Group, Trajectory, FC, SD) %>%
  xtable() %>%
  kable() %>%
  kable_styling(bootstrap_options=c("striped", "hover"))

Group	Trajectory	FC	SD
Control	Resilient	0.4696456	0.2361862
PTSD	Chronic	0.1381222	0.1583700
PTSD	Delayed	0.1253881	0.1424886
PTSD	Recovery	0.1451220	0.1522736

As expected, the the difference between simulated control groups (Traumatized controls and Traumatized PTSD) is significant to a two-sided, unequal variance t-test.

t.test(ChunkSimilarity ~ Group, 
       aconn, 
       var.equal=F) %>%
  tidy() %>%
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

estimate	estimate1	estimate2	statistic	p.value	parameter	conf.low	conf.high	method	alternative
0.3290981	0.4696456	0.1405475	138.1598	0	17514.86	0.3244291	0.3337671	Welch Two Sample t-test	two.sided

No difference between trajectories in the PTSD group

traumatized_group_aconn <- aconn %>%
  filter(Group != "Control")

traumatized_group_aconn %>%  
  aov(ChunkSimilarity ~ Trajectory, .) %>%
  summary() %>%
  xtable() %>%
  kable() %>%
  kable_styling(bootstrap_options=c("striped", "hover"))

	Df	Sum Sq	Mean Sq	F value	Pr(>F)
Trajectory	2	0.2654311	0.1327155	5.635889	0.0035837
Residuals	7002	164.8851241	0.0235483	NA	NA

And the pairwise comparisons:

pairwise.t.test(traumatized_group_aconn$ChunkSimilarity, 
                traumatized_group_aconn$Trajectory, 
                p.adjust.method = "bonferroni", paired=F, pool.sd = F) %>%
  tidy() %>%
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

group1	group2	p.value
Delayed	Chronic	0.1119543
Recovery	Chronic	0.2482275
Recovery	Delayed	0.0019504

Analysis of PTSD Model Simulations

Andrea Stocco

2/7/2021

Model and Simulations

Trajectories to simulate

Underlying memory model

Neuroiological interpretation

Event Distribution

Load and Transform the Data

Recovery Trajectories

Chronic Probability of Traumatic Intrusions

Finding a baseline to compare the model to existing literature

Factors affecting the distribution of trajectories

External factors: PTE Intensity and Context Similarity

Effects on Timecourses

Effects on Trajectories

Resilient Trajectory

Recovery Trajectory

Chronic Trajector

Delayed Trajectory

Internal Factors: Working Memory and Recollection Vividness

Effects on the Timecourse of Intrusive Memories

Effects on Trajectories

Resilient Trajectory

Chronic Trajectory

Recovery Trajectory

Delayed Trajectory

Comorbidty: Effects of Rumination

Effects on the Timecourse of Intrusive Memories

Effects on Trajectories

Resilient Trajectory

Chronic Trajectory

Recovery Trajctory

Neurobiology of PTSD

Hippocampus Size

Visualizing Hippocampus Volume Decrease

Comparisons

Functional Connectivity

Visualizing Functional Connectivity

Comparisons